#1. install necessary libraries
library(Seurat)
library(dplyr)
library(ggplot2)
library(patchwork)
options(future.globals.maxSize = 3 * 1024^3) # Set to 3 GiB
#2 Load Control 10x Data
control = Load10X_Spatial(
"/Users/cougerjaramillo/Desktop/MS/Intersession 2026/Applications/counts_and_images/1-1",
filename = "filtered_feature_bc_matrix.h5",
assay = "Spatial",
slice = "slice1",
bin.size = NULL,
filter.matrix = TRUE,
to.upper = FALSE,
image = NULL
)
# Set orig.ident of control to 0
control@meta.data$orig.ident <- "0"
SpatialFeaturePlot(control, features = "nCount_Spatial") + theme(legend.position = "right")
control
An object of class Seurat
36601 features across 3742 samples within 1 assay
Active assay: Spatial (36601 features, 0 variable features)
1 layer present: counts
1 spatial field of view present: slice1
AD = Load10X_Spatial(
"/Users/cougerjaramillo/Desktop/MS/Intersession 2026/Applications/counts_and_images/2-8",
filename = "filtered_feature_bc_matrix.h5",
assay = "Spatial",
slice = "slice2",
bin.size = NULL,
filter.matrix = TRUE,
to.upper = FALSE,
image = NULL
)
# Set orig.ident of AD to 1
AD@meta.data$orig.ident <- "1"
SpatialFeaturePlot(AD, features = "nCount_Spatial") + theme(legend.position = "right")
AD
An object of class Seurat
36601 features across 3445 samples within 1 assay
Active assay: Spatial (36601 features, 0 variable features)
1 layer present: counts
1 spatial field of view present: slice2
#3. Pre-processing
#control
#Mitochondrial DNA counts
control[["percent.mt"]] <- PercentageFeatureSet(control, pattern = "^MT-")#
#Visualize
VlnPlot(control, features = c("nFeature_Spatial", "nCount_Spatial", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(control, feature1 = "nCount_Spatial", feature2 = "percent.mt")
plot2 <- FeatureScatter(control, feature1 = "nCount_Spatial", feature2 = "nFeature_Spatial")
plot1
plot2
#Subset for high quality data
control <- subset(control, subset = nFeature_Spatial > 500 & nFeature_Spatial < 6000 & percent.mt < 40)
#Log normalize
control <- NormalizeData(control, normalization.method = "LogNormalize", scale.factor = 10000)
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
#Feature Selection
control <- FindVariableFeatures(control, selection.method = "vst", nfeatures = 2000)
Calculating gene variances
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(control), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(control)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1
plot2
#Scale Data
all.genes <- rownames(control)
control <- ScaleData(control, features = all.genes)
|
| | 0%
|
|= | 3%
|
|== | 5%
|
|=== | 8%
|
|==== | 11%
|
|===== | 14%
|
|====== | 16%
|
|======= | 19%
|
|======== | 22%
|
|========= | 24%
|
|========== | 27%
|
|=========== | 30%
|
|============ | 32%
|
|============= | 35%
|
|============== | 38%
|
|=============== | 41%
|
|================ | 43%
|
|================= | 46%
|
|================== | 49%
|
|================== | 51%
|
|=================== | 54%
|
|==================== | 57%
|
|===================== | 59%
|
|====================== | 62%
|
|======================= | 65%
|
|======================== | 68%
|
|========================= | 70%
|
|========================== | 73%
|
|=========================== | 76%
|
|============================ | 78%
|
|============================= | 81%
|
|============================== | 84%
|
|=============================== | 86%
|
|================================ | 89%
|
|================================= | 92%
|
|================================== | 95%
|
|=================================== | 97%
|
|====================================| 100%
#check final object for QC
control
An object of class Seurat
36601 features across 3657 samples within 1 assay
Active assay: Spatial (36601 features, 2000 variable features)
3 layers present: counts, data, scale.data
1 spatial field of view present: slice1
#AD
#Mitochondrial DNA counts
AD[["percent.mt"]] <- PercentageFeatureSet(AD, pattern = "^MT-")#
#Visualize
VlnPlot(AD, features = c("nFeature_Spatial", "nCount_Spatial", "percent.mt"), ncol = 3)
plot1 <- FeatureScatter(AD, feature1 = "nCount_Spatial", feature2 = "percent.mt")
plot2 <- FeatureScatter(AD, feature1 = "nCount_Spatial", feature2 = "nFeature_Spatial")
plot1
plot2
#Subset for high quality data
AD <- subset(AD, subset = nFeature_Spatial > 500 & nFeature_Spatial < 6000 & percent.mt < 40)
#Log normalize
AD <- NormalizeData(AD, normalization.method = "LogNormalize", scale.factor = 10000)
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
#Feature Selection
AD <- FindVariableFeatures(AD, selection.method = "vst", nfeatures = 2000)
Calculating gene variances
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(AD), 10)
# plot variable features with and without labels
plot1 <- VariableFeaturePlot(AD)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1
plot2
#Scale Data
all.genes <- rownames(AD)
AD <- ScaleData(AD, features = all.genes)
|
| | 0%
|
|= | 3%
|
|== | 5%
|
|=== | 8%
|
|==== | 11%
|
|===== | 14%
|
|====== | 16%
|
|======= | 19%
|
|======== | 22%
|
|========= | 24%
|
|========== | 27%
|
|=========== | 30%
|
|============ | 32%
|
|============= | 35%
|
|============== | 38%
|
|=============== | 41%
|
|================ | 43%
|
|================= | 46%
|
|================== | 49%
|
|================== | 51%
|
|=================== | 54%
|
|==================== | 57%
|
|===================== | 59%
|
|====================== | 62%
|
|======================= | 65%
|
|======================== | 68%
|
|========================= | 70%
|
|========================== | 73%
|
|=========================== | 76%
|
|============================ | 78%
|
|============================= | 81%
|
|============================== | 84%
|
|=============================== | 86%
|
|================================ | 89%
|
|================================= | 92%
|
|================================== | 95%
|
|=================================== | 97%
|
|====================================| 100%
#check final object for QC
AD
An object of class Seurat
36601 features across 3366 samples within 1 assay
Active assay: Spatial (36601 features, 2000 variable features)
3 layers present: counts, data, scale.data
1 spatial field of view present: slice2
#4. Merge the control & AD
#Merging Slices
brain.merge <- merge(control, AD, project = "brain")
#VariableFeatures(brain.merge) <- c(VariableFeatures(control), VariableFeatures(AD))
#Re-Scale Data
all.genes <- rownames(brain.merge)
brain.merge <- ScaleData(brain.merge, features = all.genes)
|
| | 0%
|
|= | 3%
|
|== | 5%
|
|=== | 8%
|
|==== | 11%
|
|===== | 14%
|
|====== | 16%
|
|======= | 19%
|
|======== | 22%
|
|========= | 24%
|
|========== | 27%
|
|=========== | 30%
|
|============ | 32%
|
|============= | 35%
|
|============== | 38%
|
|=============== | 41%
|
|================ | 43%
|
|================= | 46%
|
|================== | 49%
|
|================== | 51%
|
|=================== | 54%
|
|==================== | 57%
|
|===================== | 59%
|
|====================== | 62%
|
|======================= | 65%
|
|======================== | 68%
|
|========================= | 70%
|
|========================== | 73%
|
|=========================== | 76%
|
|============================ | 78%
|
|============================= | 81%
|
|============================== | 84%
|
|=============================== | 86%
|
|================================ | 89%
|
|================================= | 92%
|
|================================== | 95%
|
|=================================== | 97%
|
|====================================| 100%
#Visualize Results
plot1 <- VlnPlot(brain.merge, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(brain.merge, features = "nCount_Spatial") + theme(legend.position = "right")
plot1
plot2
#5 # Dimensionality reduction, clustering, and visualization
#PCA
brain.merge <- RunPCA(brain.merge, assay = "Spatial", verbose = FALSE)
Error: vector memory limit of 16.0 Gb reached, see mem.maxVSize()
#6 Identifying Cortical Layers by Heatmpap of canonical markers
# Define the list of known markers organized by layer
known_markers_list <- list(
Layer1 = c("SPARC", "MALAT1", "CXCL14", "ADIRF"),
Layer2 = c("ENC1", "HPCAL1", "CALB2"),
Layer3 = c("HOPX"),
Layer4 = c("RORB", "NUAK1"),
Layer5 = c("TUBA1B", "TMSB10", "SYT1", "STMN2", "MAP1B", "SNAP25", "PCP4", "UCHL1", "CLSTN2", "SNCG", "SMYD2", "RTN1", "STMN1", "RTN3", "NEFL", "SNRPN", "BASP1", "SYN2", "CALM3", "ENO2", "SNCA", "GAP43", "NAPB", "ELAVL4", "FXYD6", "NDRG4", "NRSN1", "RAB3C", "UHMK1", "SNAP91", "ATP6V1B2", "SLC25A22"),
Layer6 = c("MAP2K1", "DIRAS2", "KRT17"),
LayerWM = c("MBP", "MOBP", "CLDND1", "BCAS1", "MTURN", "PAQR6", "HIPK2", "DYNC1LI2")
)
# Combine all markers into a single vector and create a corresponding layer vector
combined_markers <- unlist(known_markers_list)
layer_labels <- rep(names(known_markers_list), times = sapply(known_markers_list, length))
# Ensure combined markers exist in the Seurat object
valid_markers <- combined_markers[combined_markers %in% rownames(brain.merge)]
valid_layer_labels <- layer_labels[combined_markers %in% rownames(brain.merge)]
# Check if we have valid markers before plotting
if (length(valid_markers) > 0) {
# Create a named vector for y-axis annotations
names(valid_layer_labels) <- valid_markers
# Generate the heatmap
heatmap_plot <- DoHeatmap(brain.merge, features = valid_markers) +
NoLegend() +
scale_y_discrete(labels = valid_layer_labels) + # Update y-axis labels
ggtitle("Combined Heatmap of Known Markers by Layer") # Title for the heatmap
# Display the heatmap
print(heatmap_plot)
} else {
cat("No valid markers found for heatmap generation.\n")
}
#7. Subset the Data by Region
#Subset out anatomical regions
## Subset to clusters of interest
WM <- subset(brain.merge, idents = c(3, 7))
# Extract and attach spatial coordinates from image 1
centroids <- WM[["slice1"]]@boundaries$centroids
coords <- setNames(as.data.frame(centroids@coords), c("x", "y"))
rownames(coords) <- centroids@cells
WM$x <- coords[colnames(WM), "x"]
WM$y <- coords[colnames(WM), "y"]
# After subsetting, we renormalize WM
WN <- NormalizeData(WM, normalization.method = "LogNormalize", scale.factor = 10000)
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
WM <- RunPCA(WM, verbose = FALSE)
SpatialDimPlot(WM, label = TRUE, label.size = 3)
#8 WM DEGs
# Join Layers
WM <- JoinLayers(WM)
# Set the identity classes based on the 'orig.ident' column
WM <- SetIdent(WM, value = WM$orig.ident)
# Compare cells with orig.ident values of 0 and 1
WM.AD.markers <- FindMarkers(WM, ident.1 = 0, ident.2 = 1)
| | 0 % ~calculating
|+ | 1 % ~01m 07s
|+ | 2 % ~01m 06s
|++ | 3 % ~01m 05s
|++ | 4 % ~01m 04s
|+++ | 5 % ~01m 03s
|+++ | 6 % ~01m 02s
|++++ | 7 % ~01m 01s
|++++ | 8 % ~59s
|+++++ | 9 % ~58s
|+++++ | 10% ~57s
|++++++ | 11% ~56s
|++++++ | 12% ~55s
|+++++++ | 13% ~54s
|+++++++ | 14% ~54s
|++++++++ | 15% ~53s
|++++++++ | 16% ~52s
|+++++++++ | 17% ~52s
|+++++++++ | 18% ~51s
|++++++++++ | 19% ~50s
|++++++++++ | 20% ~50s
|+++++++++++ | 21% ~49s
|+++++++++++ | 22% ~49s
|++++++++++++ | 23% ~48s
|++++++++++++ | 24% ~47s
|+++++++++++++ | 25% ~47s
|+++++++++++++ | 26% ~46s
|++++++++++++++ | 27% ~46s
|++++++++++++++ | 28% ~45s
|+++++++++++++++ | 29% ~44s
|+++++++++++++++ | 30% ~44s
|++++++++++++++++ | 31% ~44s
|++++++++++++++++ | 32% ~43s
|+++++++++++++++++ | 33% ~42s
|+++++++++++++++++ | 34% ~42s
|++++++++++++++++++ | 35% ~41s
|++++++++++++++++++ | 36% ~41s
|+++++++++++++++++++ | 37% ~40s
|+++++++++++++++++++ | 38% ~39s
|++++++++++++++++++++ | 39% ~39s
|++++++++++++++++++++ | 40% ~38s
|+++++++++++++++++++++ | 41% ~37s
|+++++++++++++++++++++ | 42% ~37s
|++++++++++++++++++++++ | 43% ~36s
|++++++++++++++++++++++ | 44% ~35s
|+++++++++++++++++++++++ | 45% ~35s
|+++++++++++++++++++++++ | 46% ~34s
|++++++++++++++++++++++++ | 47% ~33s
|++++++++++++++++++++++++ | 48% ~33s
|+++++++++++++++++++++++++ | 49% ~32s
|+++++++++++++++++++++++++ | 50% ~31s
|++++++++++++++++++++++++++ | 51% ~31s
|++++++++++++++++++++++++++ | 52% ~30s
|+++++++++++++++++++++++++++ | 53% ~30s
|+++++++++++++++++++++++++++ | 54% ~29s
|++++++++++++++++++++++++++++ | 55% ~28s
|++++++++++++++++++++++++++++ | 56% ~28s
|+++++++++++++++++++++++++++++ | 57% ~27s
|+++++++++++++++++++++++++++++ | 58% ~27s
|++++++++++++++++++++++++++++++ | 59% ~26s
|++++++++++++++++++++++++++++++ | 60% ~25s
|+++++++++++++++++++++++++++++++ | 61% ~25s
|+++++++++++++++++++++++++++++++ | 62% ~24s
|++++++++++++++++++++++++++++++++ | 63% ~23s
|++++++++++++++++++++++++++++++++ | 64% ~23s
|+++++++++++++++++++++++++++++++++ | 65% ~22s
|+++++++++++++++++++++++++++++++++ | 66% ~21s
|++++++++++++++++++++++++++++++++++ | 67% ~21s
|++++++++++++++++++++++++++++++++++ | 68% ~20s
|+++++++++++++++++++++++++++++++++++ | 69% ~19s
|+++++++++++++++++++++++++++++++++++ | 70% ~19s
|++++++++++++++++++++++++++++++++++++ | 71% ~18s
|++++++++++++++++++++++++++++++++++++ | 72% ~17s
|+++++++++++++++++++++++++++++++++++++ | 73% ~17s
|+++++++++++++++++++++++++++++++++++++ | 74% ~16s
|++++++++++++++++++++++++++++++++++++++ | 75% ~15s
|++++++++++++++++++++++++++++++++++++++ | 76% ~15s
|+++++++++++++++++++++++++++++++++++++++ | 77% ~14s
|+++++++++++++++++++++++++++++++++++++++ | 78% ~14s
|++++++++++++++++++++++++++++++++++++++++ | 79% ~13s
|++++++++++++++++++++++++++++++++++++++++ | 80% ~12s
|+++++++++++++++++++++++++++++++++++++++++ | 81% ~12s
|+++++++++++++++++++++++++++++++++++++++++ | 82% ~11s
|++++++++++++++++++++++++++++++++++++++++++ | 83% ~11s
|++++++++++++++++++++++++++++++++++++++++++ | 84% ~10s
|+++++++++++++++++++++++++++++++++++++++++++ | 85% ~09s
|+++++++++++++++++++++++++++++++++++++++++++ | 86% ~09s
|++++++++++++++++++++++++++++++++++++++++++++ | 87% ~08s
|++++++++++++++++++++++++++++++++++++++++++++ | 88% ~07s
|+++++++++++++++++++++++++++++++++++++++++++++ | 89% ~07s
|+++++++++++++++++++++++++++++++++++++++++++++ | 90% ~06s
|++++++++++++++++++++++++++++++++++++++++++++++ | 91% ~06s
|++++++++++++++++++++++++++++++++++++++++++++++ | 92% ~05s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 93% ~04s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94% ~04s
|++++++++++++++++++++++++++++++++++++++++++++++++ | 95% ~03s
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96% ~02s
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~02s
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~01s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~01s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01m 02s
# View the results
head(WM.AD.markers)